#include <iostream>
#include <fstream>
#include <string>
#include <sstream>

#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 2>::tiny floatTinyArray_t;
typedef TunaArray<double, 2>::huge ScalarField2D;

void writeTofile_DX(floatTinyArray_t *, int, string);

int main () 
{
  double length_x, length_y, dt, left_extreme, right_extreme;
  int num_nodes_x, num_nodes_y, Npoints, Nsteps, frec, Ncicles ;    

// Reading data from the input file
    std::ifstream input_cin ("input02");
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> left_extreme
	      >> right_extreme
	      >> Npoints
	      >> Nsteps
	      >> dt
	      >> Ncicles
	      >> frec;
    input_cin.close();

    StructuredMesh< Uniform<double, 2> > mesh(length_x, num_nodes_x, 
					      length_y, num_nodes_y);
    floatTinyArray_t deltas = mesh.getDeltas();

    double factor = (right_extreme - left_extreme) / (Npoints - 1);
    floatTinyArray_t line[Npoints];

    cout << "\n ----- Line construction -----\n\n";
    
    for(int i = 0; i < Npoints; i++) {
	line[i](0) = left_extreme + i * factor;
	line[i](1) = 0.5;
    }

    writeTofile_DX(line, 0, "line.");

    cout << "\n ----- Tracking: begins ----- \n\n";

// Velocity field at time t
    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());

// Velocity field at time t + dt
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());

    ParticleTracking< RungeKutta4< Linear<double, 2> > 2> particle;
    
    int count = 1;
    for(int cicles = 1; cicles <= Ncicles; ++cicles) {
	cout << " ----> Cicle " << cicles << endl;
	for(int step = 1; step <= Nsteps; ++step) {
	    cout << " step = " << step 
		 << " count = " << count << endl;

// Read the velocity at time t from a file
	    InOut::readFromFile_1(u, v, step, "velc.");

// Read the velocity at time t + dt from a file
	    InOut::readFromFile_1(un, vn, step + 1, "velc.");

// Track each particle on the line
	    for(int i = 0; i < Npoints; i++) {
		line[i] = particle.track(line[i], deltas, u, v, un, vn, dt);
	    }
	    if( !(count % frec) )
		writeTofile_DX(line, count, "line.");
	    count++;
	}
    }
    
    cout << "\n ----- Tracking: ends ----- \n\n";
    return 0;
}

void writeTofile_DX(floatTinyArray_t *line, int n, string filename)
{
    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    cout << "writing " << filename << endl;
    
    file << " object 1 class array type float rank 1 shape 2 items "
	 << 50000 << " data follows\n";
    for(int i = 0; i < 50000; ++i)
	file << line[i](0) << "   " << line[i](1) << "   ";
    file << endl;

    file << "object 2 class array type int rank 1 shape 2 items "
	 << 50000 - 1 << " data follows\n";
    for(int i = 0; i < 50000 - 1; ++i)
	file << i << "   " << i + 1 << "   ";
    file << endl;

    file << "attribute \"element type\" string \"lines\"\n"
	 << "attribute \"ref\" string \"positions\"\n";
    file << "object 3 class array type float rank 0 items "
	 << 50000 << " data follows\n";
    for(int i = 0; i < 50000; ++i)
	file << 1 << "\t";
    file << endl;
 
    file << "attribute \"dep\" string \"positions\"\n";
    file << "object \"irregular positions regular connections\" class field\n"
	 << "component \"positions\" value 1\n"
	 << "component \"connections\" value 2\n" 
	 << "component \"data\" value 3\n" 
	 << "end\n" ;
    file.close();
}


